1 Introduction

This exercise will build on the skills acquired in the Mapping Basics and Species Distribution Modeling exercises to perform Home Range Analysis and map the results. We will start by examining the structure of data used in home range analyses, we will look at ways to perform QA/QC checks on the location data, and finally various ways to visualize the data. One of the new skills we will discuss in this exercise is to create functions and loops to replicate analyses.

1.1 Packages used in this exercise

There are several specialty packages that will be used in this exercise due to the specific nature of the analyses. Some of these packages you will need to install while several others we have used in previous exercises and should already be installed.

To begin, we will install the following:

library(adehabitatHR)
library(data.table)
library(ggfortify)
library(grid)
library(move)
library(moveVis)
library(pbapply)
library(plotly)
library(rgdal)
library(sp)
library(tidyverse)
library(viridis)
library(ggplot2)
library(sf)
library(raster)
library(rasterVis)
library(leaflet)
library(dismo)
library(maps)
library(ggspatial)
library(prettydoc)
library(readxl)

2 Dataset

Let’s import a data set that tracked wandering albatross individuals with GPS devices. We will need to filter through some of it to reduce it to only two individuals since it is such a large file with many data points.

We will be looking at the Wandering Alabatross today!

A dude stands by an alabatross

Distribution of wandering albatross babies.

Import Dataset with read.csv

datas <- read.csv("./Data/alb.csv")

data <- subset(datas, id == "2016112" | id == "16-0403")
head(data)

Now we have a data set loaded in that contains information like time and date(timestamp), lat/long, and individual ID. In order to be sure that the dataset contains no outliers, we can plot the data using ggplot and plotly to interactively view the data.

qaqc_plot <- ggplot() + geom_point(data=data, 
                                   aes(lon,lat,
                                       color=id)) +
                        labs(x="Longtitude", y="Latitude") +
                        guides(color=guide_legend("Identifier"))

ggplotly(qaqc_plot)

##Reference

colors <- c("brown","yellow","darkgreen")

leaflet() %>% 
  addTiles() %>%
  addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
  addProviderTiles(providers$Esri.NatGeoWorldMap, group = "NatGeo") %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "ESRI") %>%
    addCircleMarkers(data$lon,
                   data$lat,
                   popup = data$id,
                   weight = 2,
                   color = "pink",
                   fillOpacity = 0.2,
                   group= "16-0403") %>%
  addMiniMap(position = 'topright',
             width = 100, 
             height = 100,
             toggleDisplay = FALSE) %>%
  addScaleBar(position = "bottomright")

This is just a leaflet map I added in to help show where this data is located. You will notice I am in the ocean, so the maps will suck and I apologize.

With plotly, similar to leaflet, we have the ability to examine the spread of the data points and additional information from various columns. From this plot we can see that there are two individuals, 16-0403 and 2016112.

So here we will create a function and use the lapply command to apply a function over a list or vector dataset. Specifically, this function will take the original dataset, split it into separate files based on the individual identifier, and create new *.csv files using the identifier as the filename.

lapply(split(data, data$id), 
       function(x)write.csv(x, file = paste(x$id[1],".csv", sep = ""), row.names = FALSE))

The anatomy of the function is as follows:

lapply(), apply the function over a list
split(), separates the data
function(), compose a series of steps to be applied to the data
write.csv(), write a csv file
paste(), create a character string to be used for the file name

You will now be able to see the two .csv files in the same folder of the project. Let’s create “file”, a list of the two individuals’ files.

file<-c("16-0403.csv", "2016112.csv")

3 Analysis

3.1 Imagery

As we have in previous exercises we will use raster imagery to provide additional spatial detail to the analysis. We will use map_data for this example.

In order to create a bounding box for the imagery we will use the min and max values from the coordinates in the dataset to compute those points.

world_map <- map_data("world")

ggplot() + 
           geom_polygon(data = world_map, aes(x=long, y = lat, group = group), fill = "gray", color="black") +  
  geom_point(data=data, aes(x = lon, y = lat, color = id), size = 4, alpha = 0.8) +
    coord_fixed(xlim = c(20, 80),  ylim = c(-64, -30), ratio = 1.2) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Tracking the Wandering Albatross")

3.2 Home Range Analysis

There are three basic types of home range analyses we will perform in this exercise: Minimum Convex Polygon (MCP), Kernel-Density Estimation (KDE), and Brownian Bridge Movement Model (BB). There are a number of different tuning parameters that can be applied these analyses, however in this exercise we will use the most basic versions of the analysis.

In the section above we use the lapply command to loop a function used to separate the original dataset into individual files. This is a useful tool, however, when the function loops through dozens or even hundreds of files, the process can take a long period of time to complete. Using the pblapply command adds a progress bar (i.e. pb) to the process which provides an estimated time for completion of the function. We will use a similar process to the one above, using the pblapply command to run MCP analysis on the individual *.csv files.

3.2.1 Minimum Convex Polygon

A description of the steps within the following code will be discussed following the output. The process works by establishing the function and then running pblapply referring back to the files we created in the dataset portion of this exercise.

mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$lon)
  y <- as.data.frame(data$lat)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$id)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- ggplot() + 
    geom_polygon(data = world_map, aes(x=long, y = lat, group = group), fill = "gray", color="black") +  
    coord_fixed(xlim = c(25, 74),  ylim = c(-64, -35), ratio = 1.2) +
    geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), color="blue", alpha=0.8, fill = "lightblue") +
    geom_point(data=mcp.points, aes(x=x, y=y), color="black") + 
    labs(x="Longtitude", y="Latitude", title=mcp.points$id) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(file, mcp_raster)

## [[1]]

## Warning: Use of `mcp.poly$long` is discouraged. Use `long` instead.

## Warning: Use of `mcp.poly$lat` is discouraged. Use `lat` instead.

## 
## [[2]]

## Warning: Use of `mcp.poly$long` is discouraged. Use `long` instead.

## Warning: Use of `mcp.poly$lat` is discouraged. Use `lat` instead.

The anatomy of the function above is as follows: mcp_raster <- function(filename){... for each file listed from the project directory complete the following commands

read.csv(), read the data from a given *.csv file
xy, create a coordinate file of the easting and northing data
SpatialPointsDataFrame(), project the data to UTM Zone 47
mcp(), computes home range using the Minimum Convex Polygon estimator
fortify(), turns a map into a data frame for plotting with ggplot2
grid.text(), creates annotations for the map; in this case to show area
autoplot.OpenStreetMap(), plotting for raster and vector data in ggplot2
pblapply(), referencing the files list and function

The rest of the information can be found by examing the help menu for the various commands or looking at the display options for each portion of the script.

3.2.2 Kernel-Density Estimation

The same design above can be used to create KDE plots for each individual in the dataset. Because the process is essentially the same, only portions of the function that were not previously described will be discussed.

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$lon)
  y <- as.data.frame(data$lat)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 95)
  kde.points <- cbind((data.frame(data.proj@coords)),data$id)
  colnames(kde.points) <- c("x","y","id")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- ggplot() + 
    geom_polygon(data = world_map, aes(x=long, y = lat, group = group), fill = "grey", color="black") +  
    coord_fixed(xlim = c(25, 74),  ylim = c(-64, -35), ratio = 1.2) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat),  color="black", alpha = 0.8, fill = "blue",) +
    geom_point(data=kde.points, aes(x=x, y=y)) +
    labs(x="Latitude", y="Longtitude", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(file, kde_raster)

## [[1]]

## Warning: Use of `kde.poly$long` is discouraged. Use `long` instead.

## Warning: Use of `kde.poly$lat` is discouraged. Use `lat` instead.

## 
## [[2]]

## Warning: Use of `kde.poly$long` is discouraged. Use `long` instead.

## Warning: Use of `kde.poly$lat` is discouraged. Use `lat` instead.

The anatomy of the function above is as follows: kde_raster <- function(filename){... for each file listed from the project directory complete the following commands

kernelUD(), estimation of kernel home-range
getverticeshr(), extract home-range contour

The rest of the information can be found by examing the help menu for the various commands or looking at the display options for each portion of the script.

3.2.3 Brownian Bridge Movement

Although this appears to be working backwards, in the previous two examples you have seen how to create a function that can be looped to run the analysis for a list fo files. For this analysis, we will create a conventional script for a single individual. Portions of the script that were not previously described will be discussed following the analysis.

OPHA1 <- read.csv("16.2.csv")
date <- as.POSIXct(strptime(as.character(OPHA1$timestamp),"%Y-%m-%d %H:%M:%S", tz="Asia/Bangkok"))
OPHA1$date <- date
OPHA1.reloc <- cbind.data.frame(OPHA1$lon, OPHA1$lat,
                                as.vector(OPHA1$id),
                                as.POSIXct(date))
colnames(OPHA1.reloc) <- c("x","y","id","date")
trajectory <- as.ltraj(OPHA1.reloc, date=date, id="OPHA1")
sig1 <- liker(trajectory, sig2 = 58, rangesig1 = c(0, 5), plotit = FALSE)
opha.traj <- kernelbb(trajectory, sig1 = .7908, sig2 = 58, grid = 20000)
bb_ver <- getverticeshr(opha.traj, 2)
bb_poly <- fortify(bb_ver, region = "id", 
                   proj4string = CRS("+proj=utm +zone=47+
                                     datum=WGS84 +units=m +no_defs"))
colnames(bb_poly) <- c("x","y","order","hole","piece","id","group")
bb_image <- crop(opha.traj, bb_ver, 
                 proj4string = CRS("+proj=utm +zone=47 +
                                   datum=WGS84 +units=m +no_defs"))
bb_units <- grid.text(paste(round(bb_ver$area,2)," ha"), x=0.85,  y=0.95,
                   gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
bb.plot <- ggplot() + 
    geom_polygon(data = world_map, aes(x=long, y = lat, group = group), fill = "gray", color="black") +  
    coord_fixed(xlim = c(25, 74),  ylim = c(-64, -35), ratio = 1.2) +
    geom_tile(data=bb_image, 
            aes(x=bb_image@coords[,1], y=bb_image@coords[,2],
            fill = bb_image@data$ud)) +
  geom_polygon(data=bb_poly, aes(x=x, y=y, group = group), color = "black", fill = NA) +
  scale_fill_viridis_c(option = "inferno") + annotation_custom(bb_units) +
  labs(x="Easting (m)", y="Northing (m)", title="OPHA1") +
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
bb.plot

The anatomy of the script above is as follows:

as.POSIXct(), manipulate objects to represent calendar dates and times
cbind.data.frame(), used to combine columns avoiding factorization
as.ltraj, convert data to trajectory class
liker(), used to find the maximum likelihood estimation of the parameter sig1 in kernelbb()
scale_fill_viridis_c(), create color scale for continuous data

The rest of the information can be found by examing the help menu for the various commands or looking at the display options for each portion of the script.

Try it yourself! Using the example script above, and following the example functions, create a looping function with pblapply to run the BB analysis for each individual in the files list.

3.3 Animate Trajectory Data

The final analysis we will perform in this exercise is a visual animation of the relocation data. Using the data from above you can plot(trajectory) and see a plot of the relocation information for an individual. However, we can use the move and moveVis packages to create an animation of the relocations.

We need to begin by creating a move object containing relocations (x,y), time and date information (time), a projection string, individual identifier (animal), and sensor type (sensor). See the description for move() in the help menu for optional information.

alb_order <- OPHA1[-c(14,15),]
opha.move <- move(x=alb_order$lon, 
             y=alb_order$lat, 
             time=as.POSIXct(alb_order$timestamp, 
                             format="%Y-%m-%d %H:%M:%S", tz="Asia/Bangkok"),
             proj=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"),
             data=alb_order, animal=alb_order$id)

With this information, we can now create the timing, number of frames, and animation for the relocations. First, to understand the time lag between the relocations we can use median(timeLag(opha.move, unit = "mins")) to calculate the median value. However, in the aling_move() function ther is an option to set the uniform scale manually or by using a fixed resolution: min, max, or mean. For simplicity, we will use the max temporal resolution.

movement <- align_move(opha.move, res = "max", digit = 0, unit = "mins")

With the data now on a uniform time scale we can create the frames and animation for the relocations. For this step I will use a basemap from MapBox which require token access through the use of their API. To do this you need to register with MapBox and create an access token. Then create a .Renviron file in your project folder. Copy the token information from MapBox and create an object in the .Renviron file such as map_token = 'paste token here' and add map_token = Sys.getenv('map_token') to the script below. However using the get_maptypes() script you can see there are various map services and map types that can be used. A simple output would be to use map_service = 'osm' (OpenStreetMaps) and map_type = 'topographic' or other map types available by viewing get_maptypes('owm'). when using a basemap without token access the map_token option can be removed from the script below.

Warning, the number of frames, frames/second, and output file type will determine how long this process will take and how large the output file will be.

frames <- frames_spatial(movement, path_colours = "red",
                         map_service = "mapbox", 
                         map_type = "satellite",
                         map_token = "sk.eyJ1IjoiY3VsbGV5MjAiLCJhIjoiY2t2azMwa3oyMGg2cDJvbzBrbmZsdXlxZyJ9.qHIn4uhmWFURRJuTepfATQ",
                         alpha = 0.5) %>% 
  add_labels(x = "Longitude", y = "Latitude") %>%
  add_northarrow() %>% 
  add_scalebar(distance = 2) %>% 
  add_timestamps(movement, type = "label") %>% 
  add_progress()

## Checking temporal alignment...
## Processing movement data...
## Approximated animation duration: ≈ 0.4s at 25 fps for 10 frames
## Retrieving and compositing basemap imagery...
## Assigning raster maps to frames...
## Creating frames...

animate_frames(frames, fps = 5, overwrite = TRUE,
               out_file = "./moveVis-2021fps.gif")

## Rendering animation...
## Approximated animation duration: ≈ 2s at 5 fps for 10 frames

Relocation Animation

As you can see from the result above, even a ~12sec clip with <60 frames took nearly 5min to key, obtain base imagery, assign frames, and render. So keep that in mind when creating these animations. The output from this is stored in your root directory and can be used as an animation on other programs or in html documents.

Home Range Analysis

Lexi

1 Introduction

1.1 Packages used in this exercise

2 Dataset

3 Analysis

3.1 Imagery

3.2 Home Range Analysis

3.2.1 Minimum Convex Polygon

3.2.2 Kernel-Density Estimation

3.2.3 Brownian Bridge Movement

3.3 Animate Trajectory Data